Anti-windup strategies for biomolecular control systems facilitated by model reduction theory for sequestration networks

Robust perfect adaptation, a system property whereby a variable adapts to persistent perturbations at steady state, has been recently realized in living cells using genetic integral controllers. In certain scenarios, such controllers may lead to “integral windup,” an adverse condition caused by saturating control elements, which manifests as error accumulation, poor dynamic performance, or instabilities. To mitigate this effect, we here introduce several biomolecular anti-windup topologies and link them to control-theoretic anti-windup strategies. This is achieved using a novel model reduction theory that we develop for reaction networks with fast sequestration reactions. We then show how the anti-windup topologies can be realized as reaction networks and propose intein-based genetic designs for their implementation. We validate our designs through simulations on various biological systems, including models of patients with type I diabetes and advanced biomolecular proportional-integral-derivative (PID) controllers, demonstrating their efficacy in mitigating windup effects and ensuring safety.


INTRODUCTION
Control theory has long been fundamental to the advancement of engineering systems, providing principles and methods to guide system behavior toward desired outcomes.More recently, these principles have been adopted within the growing field of synthetic biology, leading to the creation of sophisticated, biomolecular feedback controllers capable of executing complex tasks (1)(2)(3)(4).In particular, integral feedback controllers (IFCs) play a unique role among these controllers due to their ability to enable robust perfect adaptation (RPA) (5) in the system.RPA is similar to homeostasis, a recurrent principle in biology, but it is even more stringent.In essence, RPA means that a specific variable in the system is always regulated to maintain a fixed steady-state value, referred to as the setpoint, regardless of any constant disturbances or uncertainties that might be present.To this end, the antithetic integral feedback (AIF) controller was introduced to provide a chemical reaction network (CRN) that mathematically realizes integral feedback control via a sequence of chemical reactions (6).This biomolecular controller proved capable of delivering RPA not only in deterministic settings but also amid the noise inherent in stochastic settings at the population level.Furthermore, it was also shown that the AIF controller indeed represents the minimal design which is both necessary and sufficient for achieving RPA in noisy environments (7,8).Building upon the theoretical groundwork, several genetic implementations of integral and quasi-IFCs were realized in various settings including in vivo (7,9,10,11,12), in vitro (13), and optogenetically in silico (14).
Following the introduction of the AIF controller, several studies were undertaken to elucidate its dynamic performance, stability, inherent trade-offs, and tuning (15)(16)(17)(18)(19)(20).To circumvent the limitations of the AIF controller, a new generation of more advanced controllers was developed.These controllers augment the AIF motif with supplementary circuitry that realizes proportional-integral-derivative (PID) controllers (21)(22)(23)(24)(25).This advanced class of controllers has been demonstrated to enhance dynamic performance and reduce noise, embodied as cell-to-cell variability, all while preserving the RPA property.Alternative strategies have also been explored to enhance the performance the AIF controller.These include using molecular buffering (26) to enhance stability, integrating ultrasensitive modules (27,28) to amplify controller gain and reduce steady-state error when the effect of dilution (7,20) is substantial, and distributing controller components across multiple cellular populations (29)(30)(31) to improve modularity and reduce burden induced by shared resources.
Despite their promise, there are still challenges that limit the full potential of biomolecular IFCs.Among these challenges, one of the most noteworthy is integral windup (32)-a phenomenon where the integral term in the controller accumulates an error over time, causing the controller to overshoot (or undershoot) the setpoint when the actuation (or sensing) is saturating.This problem can lead to extended periods of poor performance and even system instability.It becomes especially relevant in the context of biomolecular systems where saturation is often encountered due to positivity of molecular concentrations, promoter saturation, and limited resources such as enzyme concentrations and energy molecules (33)(34)(35)(36).Integral windup is not a challenge exclusive to biomolecular systems.Its origins and subsequent mitigation strategies can be traced back to the 1930s (37) in industrial applications.Actuators, with their saturating upper and lower operational limits, have long been recognized as major contributors to integral windup manifesting as poor dynamic performance or even instability.In particular, a positivity constraint on actuators is nothing but a special case of actuator saturation where the lower bound is inherently zero.This specific form of saturation, associated with the zero bound, has been acknowledged for decades as a potential catalyst for instability arising from integral windup.Classic scenarios include valves incapable of being more than fully open or fully closed (38)(39)(40) or heaters incapable of cooling (41).This phenomenon is universal, reinforcing the pervasive relevance and importance of the integral windup issue across diverse domains.This is especially more pronounced in biology where challenges imposed by promoter saturation, rate-limiting steps, and shared limited resources (33,34,42,43) in gene regulatory networks, receptor-ligand saturation (44,45) in signaling pathways, and enzyme saturation (46) in metabolic pathways are common.For instance, IFCs show considerable promise in therapeutic applications where precise regulation of biological variables is crucial.Instabilities in such settings manifest as deviations from tolerable ranges, resulting in severe consequences.This is illustrated through a simulation study of glucose regulation in patients with diabetes, conducted using a US Food and Drug Administration (FDA)-approved mathematical model (refer to Numerical simulations) (47,48).To this end, when used in any practical application, IFCs should always be accompanied by anti-windup mechanisms (37).
For a better understanding of integral windup within the framework of biomolecular controllers, consider the visual illustration depicted in Fig. 1A.This schematic portrays a closed-loop network consisting of an arbitrary process to be regulated and a feedback controller network incorporating an IFC.In this setting, the output variable of interest is species Y, which is sensed by the controller network.Subsequently, the controller network operates on the process, aiming to robustly navigate the levels of Y toward a specific target setpoint despite constant disturbances and uncertainties in the process.In some cases, the designated setpoint may demand a degree of actuation that exceeds the capacity of the actuator, due to unavoidable saturation constraints.This leads to a situation where the controller persists in attempting to achieve an unattainable level of actuation and, as a result, certain controller species accumulates.
In mathematical terms, consider the IFC represented by a variable denoted by z ∝ ∫ t 0 e(τ)dτ , where e(t) ≜ r − y(t) is the error signal denoting the current deviation of the output from the desired setpoint r.If the control action u(t) reaches saturation, thereby inhibiting the controller's ability to drive the error to zero, then the integral continues to grow and thus leading to windup.Note that a similar scenario may also happen when the sensor saturates.Therefore, it becomes essential to design a biomolecular protection circuit implementing an anti-windup strategy.This circuit should be capable of intervening to mitigate such phenomena when they occur while remaining passive when the system is operating within a safe regime.In other words, the controller should be capable of operating in two modes: a normal mode where integral feedback facilitates RPA and a protection mode where aiming to achieve RPA is forfeited as a necessary compromise to uphold safety measures.
Figure 1B graphically demonstrates how windup could occur due to an intense transient disturbance beyond the system's tolerance.This can lead to a saturation in the control action u, resulting in an accumulation in the control variable z.Consequently, it causes a persistent error e that lingers long after the disturbance has subsided, thus taking an extended period to return to zero. Figure 1C depicts the impact of implementing an anti-windup strategy.Essentially, this strategy prevents the accumulation of the control variable z beyond a defined safe threshold by "forgetting" to integrate when operating in the danger zone.As a result, it incurs a nonzero steady-state error but expedites the return to zero error once (A) A graphical portrayal of integral windup in a biomolecular system.the closed-loop network comprises an arbitrary process in a feedback configuration with a controller network involving integral feedback control.the controller's objective is to drive the regulated output Y to a prescribed setpoint r at steady state, despite constant disturbances and uncertainties-a property referred to as RPA. the controller operates by sensing the output Y, computing the integral of the error e(t) ≜ r − y(t), and feeding the integral control action back into the process.At steady state, ̇z = 0 and thus the error converges to zero.however, saturations in the actuator and/or sensor might obstruct this process, preventing the error from reducing to zero.this results in the integral of the error accumulating, leading to windup.this causes certain molecular species, denoted as Z, to increase uncontrollably, symbolized by the spiral "winding up" due to bouncing off of a saturation wall.therefore, equipping the controller with anti-windup circuitry becomes indispensable.(B) A typical behavior in the absence of anti-windup circuits.tolerable disturbances Δ can be promptly rejected by the controller, ensuring RPA.however, when the disturbance exceeds the manageable threshold, actuation u hits its saturation limit, and a control species Z starts to accumulate, triggering windup.even when the disturbance subsides, the error takes an extended period to return to zero because z(t) is substantially beyond the safe zone, resulting in a prolonged "unwinding" process.(C) Anti-windup circuits offer a solution by stepping in only when necessary, preventing the controller from reaching a dangerous state.this intervention introduces a smaller steady-state error, helps avoid the windup process, and ensures that, once the severe disturbance subsides, the error can quickly return to zero.
the disturbance recedes.In a sense, windup could be conceptualized as a "trauma" inflicted on the system due to a severe disturbance.The ramifications of this trauma tend to linger, taking a substantial amount of time to fade even after the cause of the trauma has subsided.In this context, an anti-windup strategy can be seen as a therapeutic intervention, helping the controller to move past its "traumatic" experience by promoting a form of forgetfulness toward the disruptive event.
Various anti-windup strategies have been proposed in the literature of traditional control systems (37,49,50).These strategies modify the integrator's behavior within the controller upon detecting saturation, effectively preventing excessive error accumulation and its detrimental consequences.Nonetheless, implementing these strategies in the context of biomolecular systems poses unique challenges, primarily due to the inherent complexity and constrained structure of nonlinear CRNs.Here, we address this challenge by proposing CRN-based biomolecular realization of three anti-windup topologies.We start by phenomenologically modeling the dynamics of each topology to motivate their effectiveness in preventing windup.We then present realizations of these topologies as CRNs.We propose genetic realizations of the antiwindup circuitry-a step forward toward their implementation in living cells.

Notation
Uppercase bold letters, e.g., Y, are reserved for species names.Their corresponding lowercase letters, e.g., y(t), represent their deterministic concentrations as time-varying signals, where t is the time.Over-bars, e.g., y ≜ lim t → ∞ y(t) , represent steady-state values when they exist.A tilde, e.g., ỹ(t) ≜ y(t) − y , represents the deviation from the steadystate value, and a hat, e.g., ŷ(s) , represents the Laplace transform of ỹ(t) , where s is the Laplace variable.The s and t variables are suppressed when they are clear from context.ℝ + and ℝ − denote the set of nonnegative and nonpositive real numbers, respectively.e i is a vector, of the appropriate size, whose entries are all zeros except the ith entry being 1.Let f be a function, then range( f ) denotes its range.Let g be another function, then f ∘ g denotes the composition of the two functions, that is, ( ].If z is a scalar or a time varying signal, then z + ≜ max (z, 0) and z − ≜ max (−z, 0).‖.‖ denotes the  2 norm.

RESULTS
We begin by laying out the theoretical foundation, including the introduction of key definitions pertaining to the process to be controlled.We also provide concrete examples to clarify the definitions.We subsequently close the loop using the classical IFC suffering from possible actuator and sensor saturations and extend this discussion to the biomolecular AIF controller, drawing a rigorous mathematical link to the classical IFC.We then demonstrate how windup may arise with AIF control even in the absence of molecular saturations due to the positive nature of the AIF motif.We show how this observation fits the classical IFC framework.Further, we delve into the occurrence of windup due to molecular saturations and present three different topologies, complemented by suitable CRN realizations, that effectively implement antiwindup strategies.Potential genetic implementations are also proposed.Last, we demonstrate the efficacy of our anti-windup designs using three numerical simulations, including an FDA-approved model for glucose-insulin response in patients with diabetes (type I) (47,48).

Open-loop system
Consider a general process  Δ to be controlled, depicted in Fig. 2A, which is modeled as a nonlinear dynamical system with input and output signals u and y, respectively.The process is parameterized by Δ which can be thought of as an external disturbance or uncertainty.For simplicity, we consider processes that are single-input singleoutput throughout the paper, that is, u, y ∈ ℝ; however, this can be easily generalized to multiple-input multiple-output systems.An Process or open-loop system.the input and output of the process are respectively denoted by u and y, while Δ denotes a disturbance or uncertainty.the dynamics of the process are described by a (nonlinear) dynamic operator Δ, and its input/output steady-state map for a constant disturbance or uncertainty Δ is denoted by  Δ .(B) Gene expression as an example of a simple process.the process is composed of two molecular species X 1 and X 2 representing the mRnA and proteins whose degradation rates are denoted by γ 1 and γ 2 , respectively.note that the translation rate is denoted by k 1 , and Δ represents a disturbance entering the dynamics as a basal transcription rate which could be due to a mutation in the promoter or constitutive expression from a duplicate gene.the input to this open-loop system is an induced transcription rate denoted by u, while the output is the protein concentration y = x 2 .the dynamics and input/output steady-state map are given here.(C) closing the loop with the classical iFc.ideally, the integral controller senses the output y and computes the integral of the error e(t) ≜ r − y, which signifies the deviation of the output from the setpoint r, multiplies this integral by a gain K I , and feeds the result back to the process through the actuator.in more practical scenarios, however, sensor and actuator saturation may occur, impeding the faithful transmission of signals between the controller and the process.the equations describing the dynamics of the closed-loop system are presented in eq. 9.
example of such process  Δ is given by the following nonlinear state-space model where x ∈ ℝ n represents the internal state vector and x 0 represents the initial condition.Without loss of generality, we will proceed under the assumption that initial conditions are zero, unless specified differently.Let  Δ denote the steady-state input/output map of the process.For the nonlinear state-space realization given in Eq. 1, the steady-state input/output map is given implicitly as the following set of nonlinear algebraic equations The objective here is to design a feedback controller that endows the output y with RPA, that is, it steers y to a prescribed steady-state value or setpoint y = r , despite the presence of constant disturbances or uncertainties and regardless of the initial conditions.
Let us define four important concepts that are solely related to the controlled process and disturbance represented as  Δ .In the following definitions, we denote by and as the set of feasible inputs and disturbances, respectively.

Definition 1 (supporting input)
A supporting input, if it exists, for a given desired setpoint r and a disturbance Δ is an input value that will result in a steady-state output value equal to the setpoint r.That is, a supporting input u satisfies  Δ (u) = r and thus u ∈  −1 Δ (r).

Definition 2 (admissible setpoint)
A given setpoint r is admissible if it admits a feasible supporting input.That is, for a given setpoint r and disturbance or uncertainty Δ, ∃ , a feasible supporting input u ∈ such that  Δ (u) = r.

Definition 3 (set of admissible setpoints)
For a given process  Δ and a set of feasible inputs , the set of admissible setpoints, denoted by   Δ , , is the set of all possible admissible setpoints.That is

Definition 4 (set of admissible disturbances)
For a given admissible setpoint r, the set of admissible disturbances, denoted by  r  Δ , , , is the range of all feasible disturbances that do not destroy the admissibility of the desired setpoint r.That is Note that we drop the arguments of  and  r whenever they are clear from context.Before we provide examples, we make three remarks on the introduced definitions.

Remark 1
Admissibility is a concept that depends solely on the process to be controlled, the actuation mechanism (e.g., production, removal, and saturation) and disturbance or uncertainty; it is completely independent of the controller structure.Hence, if for a given actuation mechanism and disturbance, a desired setpoint r ∉  is not admissible, it means that there is absolutely no controller (regardless of its structure) that can deliver such a setpoint.

Remark 2
Admissibility is an algebraic steady-state property since it is ultimately related to the solvability of a set of algebraic equations.Hence, admissibility is a concept that is independent of stability, and, as such, if a setpoint is admissible, the dynamics might still be unstable.

Remark 3
Throughout the paper, we have, for the sake of simplicity, presumed that for any given desired setpoint r, the corresponding supporting input u ≜  −1 Δ (r) , if it exists, is unique.This assumption can be readily relaxed because the exact value of the supporting input is not of primary importance, as long as it is admissible and results in a stable closedloop fixed point with an output coordinate y = r.

Example 1 (gene expression)
Consider the simple model for gene expression given in Fig. 2B.The actuation and disturbance here take the form of production rates u and Δ, respectively.Both are assumed to have no saturation.Consequently, only nonnegative inputs and disturbances are feasible, that is,  =  = ℝ + .The steady-state input/output map is affine in the input u as depicted in Fig. 2B.As a result, for a given desired setpoint y = r , the supporting input is u = Requiring the supporting input to be feasible ( u ≥ 0 ) yields a condition on the setpoint given by r ≥ . Consequently, the set of admissible setpoints is given by If there is no basal transcription ( Δ = 0 ), then all nonnegative setpoints are admissible.For a given setpoint r, the set of admissible disturbances is calculated by searching for the disturbances that preserve feasibility of the supporting input (u ≥ 0) to obtain

Example 2 (stable unimolecular networks)
Consider a general unimolecular network of L species, denoted by X ≜ {X 1 , X 2 , ⋯, X L }, reacting among each other via K reaction channels.Let S Δ and λ Δ (.), respectively, denote the stoichiometry matrix and propensity function of the network in the absence of an external actuation, parameterized by a disturbance or uncertainty Δ.For unimolecular reactions, the propensity function is affine in the species concentrations, that is, Without loss of generality, let the output of interest be X L and the actuated input be X 1 where the actuation takes the form of a production reaction, that is, ∅ u(t) → X 1 and so  = ℝ + .The dynamics are thus given by Let us assume that the dynamics are Hurwitz stable, which is equivalent to assuming that the eigenvalues of S Δ W Δ have strictly negative real parts over a range of disturbances or uncertainties Δ.Then, the steady-state input/output map y =  Δ (u) is calculated by setting the time derivative in Eq. 5. to zero to obtain where α Δ is the "steady-state gain" of the network in the absence of basal expression ( b Δ = 0 ), and β Δ is the "basal offset" of the network that reflects the propagation of the basal expression rates b Δ .Note that α Δ , β Δ ≥ 0 .This arises from the fact that S Δ W Δ is a Metzler and Hurwitz matrix which means that all entries of its inverse are nonpositive (51).The steady-state input/output map is an affine function.As a result, for a given desired setpoint y = r , the supporting input is u = r − β Δ ∕α Δ .For the setpoint to be admissible, the supporting input needs to be feasible, that is, u ∈  = ℝ + ; hence, the setpoint is required to satisfy r ≥ β Δ .Consequently, the set of admissible setpoints is given by This means that for the setpoint to be admissible, it needs to be higher than the basal offset dictated by the non-actuated process.For a given setpoint r, the set of admissible disturbances is Example 1 is a special case where α Δ = k 1 ∕γ 1 γ 2 and β Δ = k 1 Δ∕γ 1 γ 2 .

Closing the loop with the classical integral controller
According to the internal model principle (52), RPA can be achieved using IFCs.To this end, consider the closed-loop system depicted in Fig. 2C, where the controlled process is now in a feedback interconnection with an integral controller.Ideally, the actuator should faithfully transmit the controller's output to the controlled process (i.e., u = v), and the sensor should faithfully transmit the output to the summation junction (i.e., w = y) to compute the deviation (or error) e of the sensed output w from the desired setpoint r.However, practically, the actuator and sensor may saturate beyond certain thresholds.Even worse, they may deform the actuation and sensed signals.Sensor and actuator saturations are modeled by the saturation blocks depicted in Fig. 2C.Integral feedback control is achieved by the integrator module which computes the time integral of the error to yield the feedback signal v(t) to the actuator.The full dynamics of the closed-loop system are given by where [u min , u max ] and [y min , y max ] represent the actuator and sensor ranges, respectively, and sat(.) is the saturation function defined, for any z, a, b ∈ ℝ with a < b, as

Ideal setting
In the ideal setting, we have the following assumptions.

Assumption 2 (no saturations)
The sensor and actuator do not exhibit any saturations, i.e., u min = y min = − ∞ and u max = y max = ∞.
Under Assumptions 1 and 2, we have and, as a result, the dynamics given in Eq. 9 boils down to This guarantees that the output y will converge at steady state to the setpoint r, independent of Δ , as long as the closed-loop system remains stable, that is, RPA is achieved.

Non-ideal setting
The fixed point of the dynamics of the non-ideal setting can be computed by setting the time derivatives in Eq. 9 to zero to obtain the following set of algebraic equations The following lemma characterizes the existence of the fixed point.

Lemma 1
Consider the closed-loop system depicted in Fig. 2C, where the actuator and sensor ranges are given as [u min , u max ] and [y min , y max ], respectively.For a given desired setpoint r and a steady-state disturbance or uncertainty Δ , the closed-loop fixed point exists with a feasible supporting input u ∈ if and only if • The setpoint is admissible, i.e., • The sensor does not saturate at steady state, i.e., r ∈ [y min , y max ].
• The actuator does not saturate at steady state, i.e.,  −1 Δ (r) ∈ u min , u max .Furthermore, the fixed point is given by The proof can be found in section S3A.Lemma 1 essentially says that if the desired setpoint is inadmissible by the process, or if either the sensor or actuator saturates, then Controlled process y =  Δ (u) e.g., ẋ = f Δ (x,u) the integral controller fails and the closed-loop dynamics exhibit no equilibrium.The only way that integral control fails in achieving RPA is by losing stability of the closed loop whether this instability is due to nonexistence of a fixed point or the existence of an unstable or unreachable fixed point.
To use Lemma 1 in a biological setting, we need to extend it into the realm of CRNs.To this end, we present in Box 1 and Fig. 3 two limit theorems that investigate the dynamics and fixed-point behaviors of complex CRNs involving fast sequestration reactions.The theorems presented offer a useful model reduction tool that is valid in both deterministic and stochastic settings.These results will be instrumental throughout this study in bridging CRN motifs with established control-theoretic architectures.

AIF control
Owing to the importance of RPA in biology, the AIF controller was introduced in (6) to realize integral control as a CRN which endows the system with the desirable RPA properties.In this section, we establish a connection between the AIF controller and the classical IFC of Fig. 2C.
The basic AIF controller is depicted in the closed-loop network of Fig. 4A, where the objective is to endow the regulated output X L with RPA.The closed-loop network is constituted of the AIF controller, composed of two species Z 1 and Z 2 , connected in a feedback configuration with an arbitrary network or process, composed of L species X 1 , X 2 ,⋯, X L .For generality, we adopt the biomolecular control paradigm introduced in ( 22), where we assume that the controller can only interact with the process species X 1 and X L for actuation and sensing, respectively.The controller network involves four reaction channels listed in Table 1.
The closed-loop dynamics are thus given by Note that since the actuation by the controller is carried out via a production reaction only which is assumed to have no saturation limit, the set of feasible inputs is  = ℝ + .The integral action can be seen by looking at the dynamics of z ≜ z 1 − z 2 given by where r ≜ μ/θ denotes the setpoint.As a result, as long as closed-loop stability is maintained, the regulated output converges to the setpoint y = x L = r regardless of the initial conditions and process uncertainties/disturbances Δ.As such, the AIF controller endows the regulated output species X L with RPA.In practice, Z 1 and Z 2 undergo dilution, which results in a leaky integral controller and imperfect adaptation.However, studies have shown that when the dilution rate is small relative to other network rates, its impact is negligible and the integral controller's performance remains effective (7,20).
Box 1.A model reduction result for CRNs with a fast sequestration reaction considers the two dynamical systems,  η and , depicted pictorially in the Fig. 3. the deterministic dynamics of each of the two systems are described by the following ordinary differential equations We make the following assumption.Assumption 3. the three functions W 1 , W 2 : ℝ L+2 + → ℝ , and F: ℝ L+2 + → ℝ L are assumed to be globally lipschitz on their domains, and their form is such that the solution of  η lies in ℝ L+2 + .the following theorem states the convergence of the dynamics, as η → ∞, of the full system  η to those of the reduced system  over any compact time interval.Theorem 1.Under Assumption 3 and as η → ∞,  η converges to  in the following senses: ∀t > 0 we have A numerical demonstration of theorem 1 is provided in fig.S7.Remark 4. in the limit as η → ∞, the concentration of at least one of the two species Z 1 or Z 2 is zero at any time.hence, the dynamics may enter a low-copy regime where stochastic fluctuations become more relevant.to this end, we prove in section S2 that theorem 1 is also valid in the stochastic setting.We also provide numerical validations in fig.S8. the following theorem, on the other hand, states the convergence of the fixed points.Theorem 2. Under Assumption 3 and as η → ∞, the limit point x, z 1 , z 2 of any convergent subsequence of nonnegative fixed points of system  η can be transformed to a fixed point of system  given by x, z 1 − z 2 .the following corollary characterizes the relationship between the existence of the fixed points of systems  η , as η → ∞ and .Corollary 1.Under Assumption 3, if system  admits a unique fixed point (x, z) , then all convergent subsequences of nonnegative fixed points of system  η have a single limit point as η → ∞ given by x, z + , z − .however, if system  admits no fixed points, then there is no convergent subsequence of nonnegative fixed points of system  η .the proofs of theorems 1 and 2 and corollary 1 can be found in section S1.
To reveal the connection between the AIF controller and the classical IFC depicted in Fig. 2C, we examine the case where the sequestration rate is large enough.In the limit as η → ∞, applying Theorem 1 [with W 1 (x, z 1 , z 2 ) ≜ μ and W 2 (x, z 1 , z 2 ) ≜ θx L ] yields a reduced model for the controller given by Introducing the intermediate variables v ≜ kz, w ≜ y, and e ≜ r − w and the integral gain K I ≜ kθ allows us to draw the block diagram of the closed-loop system depicted in Fig. 4B.This block diagram allows us to directly compare the AIF controller with the classical IFC depicted in Fig. 2C.The AIF controller is essentially the same as the classical IFC where the sensor has no saturation limits, that is, [y min , y max ] = [−∞, +∞], while the actuator saturates only from below at zero, that is, [u min , u max ] = [0, +∞].This is intuitive since the actuation is carried out via the production of X 1 , and production cannot be negative.As a result, if the desired setpoint r requires a supporting input u < 0 , then the closed-loop fixed point with the reduced controller  ∞ will cease to exist.However, for the full model with large η, the closed-loop fixed point may not cease to exist but instead suffer from a negative component z 1 < 0 (see Corollary 1) rendering it unreachable (recall that the positive orthant is invariant under the dynamics of CRNs), and as a result, the dynamics become unstable.One can think of the AIF controller  η with finite η as an integral controller coupled with "dynamic saturation" from below at zero, as compared to the case with infinite η where the saturation max(z,0) becomes static.
Note that if one has the luxury of changing the actuation mechanism or even adding more actuation mechanisms, then the issue of nonexistent or unreachable fixed point can be circumvented.For instance, by actuating via both production and degradation of X 1 , the control action in Eq. 14 becomes u = kz 1 − F(z 2 , x L )x 1 , and the set of feasible inputs and actuator range become  = ℝ .Notably, for F(z 2 , x L ) = γz 2 or F(z 2 , x L ) = γx L , we obtain (filtered) proportional-integral controllers [see (22)].
We close this section, by providing a numerical demonstration depicted in Fig. 4C where the gene expression model presented in Fig. 2B is controlled by the AIF controller.The plot shows that RPA is achieved as long as the disturbance Δ ∈ r is admissible.However, when Δ transitions to an inadmissible level, although transiently, integral windup manifests.In this scenario, z 1 attempts to go negative, but of course cannot, due to the positivity of the system.Consequently, z 2 accumulates and remains high for an extended period, even after the disturbance reverts to an admissible level.This results in poor dynamic performance induced by integral windup.

AIF control with actuation/sensing saturation
Consider now the AIF controller where catalytic production reactions may saturate as depicted in the closed-loop network of Fig. 5A.The controller network here involves the same four reaction channels described in Table 1, but the propensity functions of the actuation and sensing reactions are now replaced with kh a (z 1 ) and θh s (x L ), respectively, where h a and h s are nonlinear functions that may introduce saturation.Examples of such functions are Michaelis-Menten or Hill-type functions.The closedloop dynamics are thus given by Process

System System
Fig. 3. Model reduction for CRNs with fast sequestration reactions.the left schematic illustrates the full system  η , where η denotes a sequestration rate. in contrast, the right schematic portrays the reduced system  which accurately reflects the dynamics of  η as η → ∞. the system  η consists of two species, Z 1 and Z 2 .their production and removal rates are compactly represented by the functions W 1 and W 2 , respectively.the rest of the network (green box) is arbitrary and remains unchanged in the reduced model.essentially, when the sequestration rate is fast, the two species operate as a subtraction mechanism, given by W 1 − W 2 which is subsequently integrated in time to produce Z = ∫ (W 1 − W 2 ). the positive and negative portions of Z are isolated and relayed to the rest of the network in place of Z 1 and Z 2 , respectively.Refer to Box 1 for the complete mathematical description.The integral action can still be seen by looking at the dynamics of z ≜ z 1 − z 2 given by where, r in ≜ μ/θ.Observe that, assuming closed-loop stability, h s y = r in .As a result, the sensing function h s determines the relationship between the input and output setpoints denoted by r in and r out , respectively.Note that in the ideal setting of Fig. 4, h s is the identity function, and thus, the input and output setpoints match, that is, r in = r out = r.However, in general, the input and output setpoints satisfy the following algebraic equation or r out = h −1 s r in when h s is invertible.To this end, we refine the concept of admissible setpoints into two distinct categories: admissible input setpoints, denoted as  in , and admissible output setpoints, denoted as  out .These can be mapped to each other using Eq.19 for monotonically increasing functions h s , i.e.,  in ≜ h s ( out ).To reveal the connection with the classical IFC depicted in Fig. 2C, we examine, once again the case where the sequestration rate is large enough.In the limit as η → ∞, applying Theorem 1 [with W 1 (x, z 1 , z 2 ) ≜ μ and W 2 (x, z 1 , z 2 ) ≜ θh s (x L )] yields a reduced model for the controller given by Introducing the intermediate variables v ≜ kz, w ≜ h s (y), and e ≜ r − w and the integral gain K I ≜ kθ allows us to draw the block diagram of the closed-loop system depicted in Fig. 5B.A direct comparison with the classical IFC, depicted in Fig. 2C, reveals that the difference here is embodied in the actuator and sensor modules that are now replaced with the nonlinear functions ψ a and ψ s define in Fig. 5A.These functions not only act as saturation components but also deform the signals.The following lemma slightly generalizes Lemma 1 to the case where the saturation blocks of Fig. 2C are replaced with the monotonically increasing functions ψ a and ψ s shown in Fig. 5B.

Lemma 2
Consider the closed-loop system depicted in the block diagram of Fig. 5B, where ψ a and ψ s are strictly monotonically increasing functions that may saturate.For a given desired output setpoint r out and a steady-state disturbance or uncertainty Δ , the closed-loop fixed point exists with a feasible supporting input u ∈ if and only if • The output setpoint is admissible, i.e., • The sensor does not saturate at steady state, i.e., r in ∈ range(ψ s ).
• The actuator does not saturate at steady state, i.e.,  Furthermore, the fixed point is given by The proof can be found in section S3B.Therefore, as per Corollary 1, for the full model with large η, the nonnegative closed-loop fixed point will also not exist if any of the conditions of Lemma 2 are violated.However, the closed-loop fixed point may suffer from a negative component rendering it unreachable, and as a result, the dynamics become unstable.
Two numerical simulations are depicted in Fig. 5 (C and D) where the gene expression model presented in Fig. 2B is controlled by the AIF controller that involves actuator and sensor saturation with Figure 5C demonstrates setpoint tracking, provided the desired input setpoint, tuned by μ, is admissible (i.e., r in ∈  in ).Similarly, Fig. 5D demonstrates disturbance rejection when the disturbance is admissible (i.e., Δ ∈  r ).The sets of admissible input/output setpoints and disturbances,  in / out and  rout , respectively, are calculated in Example 3, by modifying Example 1 to factor in the effects of saturation.

Example 3 (gene expression with saturation)
Consider the simple model for gene expression given in Example 1; however, replace the actuation with u = k z 1 ∕ κ a 1 + z 1 ∕ κ a .As a result, actuation saturation constrains the feasible inputs to the range = 0, k .The steady-state input/output map remains unchanged as given in Fig. 2B, resulting in an identical supporting input: u = γ 1 γ 2 k 1 r out − Δ .To ensure the supporting input's feasibility ( 0 ≤ u ≤ k ), the set of admissible output setpoints is thus given by Subsequently, the set of admissible input setpoints r in ∈  in can be mapped from  out using Eq. 19, i.e.,  in = h s ( out ). Recalling that r in ≜ μ/θ, we can outline the gray shaded area in Fig. 5C which is given by Last, for a given output setpoint r out , the set of admissible disturbances can be computed by determining the disturbances that preserve the feasibility of the supporting input (u ∈ ) , resulting in Intuitively, when the (input/output) setpoint exceeds the admissible range, species Z 1 accumulates in a futile attempt to increase actuation u, which is already saturated.Conversely, when the output Table 1.Reactions describing the basic AIF controller.

Setpoint Reaction
setpoint falls below the admissible range, species Z 1 essentially tries to become negative to reduce actuation u and, as a result, lingers at zero and thus inflicting a buildup of Z 2 .These observations apply to both Fig. 5 (C and D).

Anti-windup strategies
In this section, we introduce various anti-windup strategies designed to mitigate the unwanted effects of windup.We begin by describing these designs from a phenomenological perspective.We then present their CRN realizations.We propose genetic implementations of these anti-windup schemes.

Anti-windup topologies
The objective here is to present three distinct strategies to alleviate the effects of windup by augmenting the basic AIF motif with extra anti-windup circuitry.These strategies are illustrated as three topologies in Fig. 6A.Although the three topologies are mechanistically different, they share the same concept: The anti-windup circuitry is only activated to halt further growth when either of the controller species, Z 1 or Z 2 , increases excessively.Note that although the three topologies perform comparably in computational terms, the particular choice of topologies is application specific and depends on practical considerations of the available genetic parts.The anti-windup circuitry for all three topologies are symmetric across Z 1 and Z 2 (blue and red switches in Fig. 6A); thus, we will only explain the anti-windup operation on Z 1 for brevity.In the case of topology I, anti-windup is accomplished by having Z 2 produce Z 1 whenever its concentration reaches a high level.As a result, the produced Z 1 species sequester the overabundant Z 2 species and thus effectively prevent it from growing any further.Mathematically, this is represented by the monotonically increasing function h 1 (z 2 ), denoting an additional production rate of Z 1 , which exhibits a switchlike behavior (depicted by the red switch in Fig. 6A, topology I).
In situations where Z 2 concentration remains low, the switch is off, allowing the circuit to function purely as an integral controller.Conversely, when Z 2 concentration exceeds a specified threshold, the switch is activated, transitioning the circuit into protection mode.While this mode relinquishes the pure integral control function, it safeguards the system from the detrimental effects of excessive growth.
An example of such a function, h 1 , is given in the gray box of Fig. 6A.Here, the production rate remains zero for z 2 ≤ β 1 .However, once the threshold β 1 is exceeded, the rate becomes linear in z 2 with a positive slope α 1 .
In the case of topology II, anti-windup is accomplished by having Z 1 degrade when its concentration reaches a high level.This strategy is represented by the monotonically increasing function h 1 (z 1 ), denoting a degradation rate of Z 1 , which exhibits a switch-like behavior (depicted by the blue switch in Fig. 6A, topology II).
In topology III, the anti-windup mechanism is achieved by having Z 1 repressing its own production once its concentration surpasses a threshold.Mathematically, this scenario is depicted by the

A B C
Example for topologies I, II: Topology III: Reference/sensor conditioning  in all topologies, the controller operates normally as an AiF controller ensuring RPA, as long as Z 1 and Z 2 levels stay below a threshold.Otherwise, the anti-windup mechanisms (red and blue switches) activate, shifting the controller to protection mode.topology i increases production of Z 1 in response to high Z 2 and vice versa.topology ii degrades Z 1 and Z 2 when their levels are high.topology iii reduces production rates of Z 1 and Z 2 when their concentrations exceed a threshold.the switches are represented by monotonic functions h 1 and h 2 .For topologies i and ii, these functions are zero until the threshold and then increase.For topology iii, they are one until the threshold and then decrease.Mathematical examples of these functions, realizable as cRns, are provided in the lower gray box (see Fig. 7).introducing the intermediate variable v, input setpoint r in , integral gain K I , and functions ψ a and ψ s , as shown in the top gray box, generates the block diagram in the bottom left, depicting the reduced dynamics similar to Fig. 5B with added anti-windup dynamics.(B) compact block diagram for topologies i and ii.these topologies use a conditional integration approach, where the controller integrates when v is within the safe range [−kβ i , kβ j ] and transitions to a leaky integrator outside this range to counter windup.topology i is associated with (i, j) = (1,2), while topology ii with (2,1).(C) Block diagram for topology iii, implementing reference/sensor conditioning.this reduces the input reference or sensor signal when v is outside the safe range [−kβ 2 , kβ 1 ].Specifically, high positive v scales down the input setpoint r in , while high negative v scales down the sensed signal w.
monotonically decreasing function h 1 (z 1 ), which attenuates the maximum production rate μ whenever z 1 exceeds the threshold.An example of such a function, h 1 , is given in the gray box of Fig. 6A.Here, the production rate persists at μ for z 1 < β 1 .However, upon crossing the threshold β 1 , the rate decays with Z 1 .Note that this topology encompasses the anti-windup circuit introduced in (53).
To connect these anti-windup strategies with classical controltheoretic methodologies, we once again look at the asymptotic limit of a large sequestration rate η.By invoking the same intermediate variables as before (repeated in the gray box of Fig. 6A for convenience), we can construct the block diagrams of the controllers illustrated in Fig. 6 (B and C).The derivations once again rely on Theorem 1, and the details can be found in section S4.Both topologies I and II embody the same control-theoretic principle similar to "conditional integration"; whereas topology III applies a concept similar to "reference/sensor conditioning" (37,50,54).
In the case of conditional integration shown in Fig. 6B, the controller operates as a pure integral controller, with r in being the input setpoint, as long as v ≜ k(z 1 − z 2 ) ∈ [−kβ i , kβ j ].However, as soon as v departs from this range, the controller behaves like a "leaky integrator, " represented by a transfer function of the form of K I s + α .This essentially incorporates a forgetting factor that "forgets" the need to integrate the entire history of the error signal.While leaky integrators introduce a steady-state error, they reciprocate with a stabilizing effect (16).The steady-state error can be shown to be proportional to the deviation of v from the integrator regime [−kβ i , kβ j ].To see this, we resort to the dynamics of v described in the block diagram of Fig. 6B that can be explicitly written as Recalling that K I ≜ kθ, the error at steady state, when it exists, can be written as and hence the steady-state error is indeed proportional to the deviation of v from the integrator regime.
In the case of reference/sensor conditioning depicted in Fig. 6C, the controller operates as a pure integral controller, with r in being the input setpoint, as long as v ≜ k(z 1 − z 2 ) ∈ [−kβ 2 , kβ 1 ].However, as soon as v ventures outside this range, either the reference or sensor signal undergoes a dynamic modification or "conditioning" based on v and thus steering a conditioned version of the error toward zero.This strategy essentially gives up trying to track the desired setpoint when necessary to prevent windup.More specifically, if v > kβ 1 , then the input reference r in is conditioned by reducing it by a factor of kα 1 /(kα 1 + v − kβ 1 ).Hence, this conditioning becomes more pronounced the further v exceeds the integrator operating regime [−kβ 2 , kβ 1 ].Similarly, if v < kβ 2 , then the sensor signal w is conditioned by reducing it by a factor of kα 2 /(kα 2 − v − kβ 2 ).As a result, the conditioning is more substantial the further v descends from the integrator operating regime.

Sequestration-based switches
To realize the three anti-windup strategies using CRNs, we first need to realize the switches depicted in Fig. 6A, which are mathematically represented as the monotonic functions h 1 and h 2 .Specifically, we consider the functional forms shown in the gray box of Fig. 6A and reproduced in Fig. 7 for convenience.The sequestration networks depicted in Fig. 7 (A and B) deliver a steady-state response that approximates the desired functions.This approximation becomes exact in the limit of strong sequestration.
The two sequestration networks consist of two species V 1 and V 2 which sequester each other at a rate of η while separately degrading at rates δ 1 and δ 2 , respectively.Here, the input, represented as u, signifies the production rate of V 2 , while the thresholds are determined by the production rate u 0 of V 1 .The distinguishing feature between the two networks lies in the output: In Fig. 7A, V 2 serves as an activator; conversely, in Fig. 7B, V 2 functions as a repressor.Accordingly, the output is represented by the activation and repression rates, respectively.Once again, Theorems 1 and 2 are invoked here, and the derivations can be found in section S5.These designs are also numerically verified in the plots shown at the bottom of Fig. 7. Notably, decreasing the degradation rates yields a steeper slope leading to an ultrasensitive response (27,55)-a feature that is not required here for anti-windup.

CRN realizations of anti-windup topologies
Next, we leverage the sequestration-based switches to fully embody the three anti-windup topologies depicted in Fig. 6A as CRNs. Figure 8 (A to C) each represents a CRN corresponding to anti-windup topologies I, II, and III, respectively.Given the symmetry of the anti-windup schemes in Z 1 and Z 2 , we describe the anti-windup operation on Z 1 only.
For the CRN realization of topology I, depicted in Fig. 8A, the red switch is implemented using the sequestration network from Fig. 7A.This network consists of two species, V 1 and V 2 , which sequester each other at a rate of η v and degrade separately at rates δ v1 and δ v2 .The production rate of V 2 , represented as g 1 (z 2 ), serves as the input to the switch and is a monotonically increasing function of z 2 .Its output, denoted by h 1 (v 2 ), is a monotonically increasing function of the activator V 2 , which plays a role in the production rate of Z 1 .Therefore, with reference to Fig. 7A and section S5, assuming a high sequestration rate η v , h 1 v 2 can be written as where v 0 is the constant production rate of V 1 , serving as the threshold.Note that h 1 and g 1 are retained as arbitrary monotonically increasing functions for generality; however, linear representations of these functions would align exactly with the depiction in Fig. 6A, as shown in the example of Fig. 8B.
For the CRN realization of topology II, the blue switch is implemented using a similar sequestration network.The primary differences lie in the production rate of V 2 , which is now a function of z 1 rather than z 2 , and V 2 now catalyzes the degradation of Z 1 at a rate h 1 (v 2 ).Thus, the degradation propensity of Z 1 is represented by h 1 (v 2 )z 1 .With reference to Fig. 7A and section S5, assuming a high sequestration rate η v , the degradation propensity at steady state can be expressed as where the propensity now grows quadratically beyond the threshold v 0 (see fig.S1B).Last, for the CRN realization of topology III, the blue switch is realized using the sequestration network from Fig. 7B.The main distinction from topology II lies in V 2 , which now represses Z 1 at a rate h 1 (v 2 ).This rate is a monotonically decreasing function of v 2 , for instance, h(v 2 ) = α 1 /(1 + v 2 /κ 1 ).

Genetic realizations of anti-windup topologies
Genetic realizations of biomolecular integral controllers have already been successfully built and tested in both bacterial and mammalian cells (1).Therefore, we will not delve into the genetic realization of integral controllers here.Instead, our focus will be on proposing genetic realizations of anti-windup circuitry.The genetic implementations proposed in this study rely on inteins.An intein is a short protein segment that can autocatalytically excise itself from a protein structure and simultaneously rejoin the remaining segments, known as exteins (56).Split inteins, a subcategory of inteins, are divided into two parts commonly referred to as IntN and IntC.When active, these split inteins can heterodimerize and independently carry out protein splicing reactions, which involve the irreversible creation and destruction of peptide bonds in a strict one-to-one stoichiometric ratio.Leveraging their ability to exchange, cleave, or ligate amino acid sequences, split inteins provide the necessary foundation to realize the sequestration reactions that are crucial to the Input/output static map Fig. 7. Biomolecular switching using sequestration motifs.the two sequestration networks in (A) and (B) consist of two species V 1 and V 2 , which sequester each other at a rate η while degrading separately at rates δ 1 and δ 2 .the input u corresponds to the production rate of V 2 , and the thresholds are determined by the production rate u 0 of V 1 .the two networks differ in the role of V 2 which acts as an activator in (A) and a repressor in (B). the bottom plots numerically confirm these designs.
) and (W 1 , W 2 ) that sequester each other, are used to implement switch-like behaviors that intervene to prevent windup when necessary.topology i uses activator switches and so h 1 and h 2 are monotonically increasing beyond the threshold.topology ii uses a similar mechanism, but the switches catalyze the selfdegradations of Z 1 and Z 2 .here, h 1 (v 2 ) and h 2 (w 2 ) represent degradation rates that are monotonically increasing beyond the thresholds.the corresponding propensity functions, however, are given by h 1 (v 2 )z 1 and h 2 (w 2 )z 2 , respectively.As a result, if h 1 and h 2 are linear functions beyond the threshold, then the corresponding propensities are quadratic.last, topology iii introduces repressor switches, and so, h 1 (v 2 ) and h 2 (w 2 ) monotonically decrease beyond the thresholds.note that while the functions h 1 , h 2 , g 1 , and g 2 are kept intentionally general in this depiction, specific examples are given in (B) and (c). in (B), the functions are presented in their linear form, while in panel (c), h 1 takes on a hill form.
genetic realization of the anti-windup CRNs depicted in Fig. 8.The flexible and small-sized nature of split inteins (57), their existence in numerous orthogonal pairs (58), and their applicability across diverse life forms (59, 60) make them particularly attractive for constructing genetic circuits.They have previously been used to build biomolecular integral controllers (9). Figure 9 (A to C) depicts three genetic realizations compatible with topologies I, II, and III from Fig. 8, respectively.Given that the anti-windup circuitry is symmetric with respect to Z 1 and Z 2 , our explanation here will focus on the genetic realization for one of them only.As depicted in Fig. 9A, the circuit is composed of two genes that express the proteins V 1 and V 2 .The first gene is driven by a constitutive promoter and includes a degradation tag fused to IntN.The second gene, on the other hand, is driven by an inducible promoter activated by Z 2 .It consists of a DNA binding domain linked via an IntC and is fused to an activation domain and a degradation tag.Upon expression of the two genes, the corresponding proteins V 1 and V 2 engage in an intein-splicing reaction and are independently degraded by the proteasome.It is important to note that before splicing, V 2 functions as an activator, driving the expression of Z 2 .However, its function is dismantled when it undergoes the splicing reaction with V 1 .The threshold of this anti-windup circuit can be adjusted by modifying the strength of the constitutive promoter or its plasmid copy number.It can also be tuned by inducers if the constitutive promoter is replaced by an inducible promoter.In contrast, the anti-windup circuit in Fig. 9C is the same as that in Fig. 9A but with only two differences: Gene 2 lacks an activation domain and its promoter is activated by Z 1 instead of Z 2 .Consequently, V 2 acts as a repressor rather than an activator, making this circuit configuration appropriate for implementing topology III.Last, the anti-windup circuit in Fig. 9B is similar to that in Fig. 9C, except for a single distinction: The DNA binding domain is replaced by a protease.This modification causes V 2 to act as a protease that degrades Z 1 rather than repressing its production.

Numerical simulations
To verify the effectiveness of our proposed anti-windup circuits, we undertake three numerical simulations using three processes of increasing complexity.Figure 10 depicts the simulations results for topology I only; however, similar results for topologies II and III can be found in figs.S4 and S5, respectively.
The first simulation, shown in Fig. 10A, revisits the gene expression process outlined in Example 3-a process that previously demonstrated windup effects as shown in Fig. 5 (C and D).Here, we use the CRN realization of topology I, shown in Fig. 8A, to mitigate these windup effects.As depicted in Fig. 10A, our anti-windup strategy successfully eliminates windup effects in both setpoint tracking (left) and disturbance rejection (right) scenarios.The inset figures track the dynamics of the anti-windup circuits V 2 and W 2 , revealing that they remain dormant until the preset thresholds of Z 1 and Z 2 are surpassed.Upon reaching these thresholds, they become active and shift the controller into a protection mode that sacrifices RPA (which is anyway unattainable) for safety, thereby preventing Z 1 and Z 2 from reaching excessive levels that would cause windup and deteriorate dynamic performance.
Next, we illustrate that our proposed anti-windup topologies are compatible with biomolecular PID controllers.In this scenario, we consider a more complex process controlled by a third-order PID controller (22), as represented in the closed-loop network of Fig. 10B.This process consists of six species, X 1 through X 6 , with X 1 being the input species and X 6 as the regulated output species of interest.Each species X i degrades at a rate γ i and catalytically produces X i+1 for i = 1,2, ⋯ ,5.In addition, X 6 actively degrades X 2 .This process is taken from (22,23) but with two modifications: Catalytic production reactions are changed to potentially saturating Hill functions rather than linear functions, and a basal expression rate b for all the species is added.The simulation results presented in Fig. 10B involve two different disturbances: a persistent but admissible disturbance in γ 6 , which the controller successfully rejects, and a transient but inadmissible disturbance in k p which is the maximum production rate of all species.The first disturbance could represent a mutation in the relevant gene's coding sequence or the activation of another pathway that sequesters the output protein.The second disturbance could model a severe, transient metabolic burden that suppresses the production of all species within the regulated network.This intense, albeit brief, disturbance triggers windup as Z 1 builds up to high levels, leading to poor dynamic performance even after the disturbance has passed.However, with the implementation of the anti-windup circuitry, the accumulation of Z 1 is circumvented, resulting in superior dynamic performance after.Once again, the inset figure displays the dynamics of the relevant anti-windup species W 2 , which only activates during windup, transitioning the controller into a protection mode.
In the final simulation example, we use a highly complex process pertaining to a mathematical whole-body model of patients with type I diabetes (47,48,61).This model gave rise to the first computer simulator approved by the FDA as an alternative to preclinical trials and animal testing.Note that a prior simulation study demonstrated an effective control of this model using the AIF controller (10), where insulin production is the control action, and blood glucose level is the output required to exhibit RPA.We retain the same AIF closed-loop network and parameter values from (10), with one exception: We increase the parameters μ and θ by a factor of 100 to speed up the response while maintaining the same setpoint at 100 mg/dl-a value within the healthy range of [70,126] mg/dl.To test for windup effects, we introduce a severe but transient disturbance in the model.Specifically, the parameter k p1 is reduced by 50% for a duration of 6 hours, simulating a marked change in endogenous glucose production (EGP), after which the EGP rate returns to normal.The gray curves in Fig. 10C correspond to the response of the model of a nondiabetic person not controlled by our circuits.These curves show that insulin secretion abruptly drops to zero to counteract the introduced severe disturbance and successfully returns the glucose level to the healthy range.It is worth noting that even in a model of a nondiabetic individual, plasma glucose transiently drops to a substantially low level, emphasizing the intensity of the applied disturbance.When we use the AIF controller without the anti-windup circuitry to control the model of a patient with diabetes, we get the response depicted in red.The accumulation of Z 1 to high levels in response to the inadmissible disturbance is followed by a lengthy recovery time-a clear sign of windup.This leads to a prolonged period of low insulin secretion as the actuator species Z 2 remains near zero, resulting in dangerously elevated glucose levels for around 20 hours.However, the implementation of our anti-windup topology I effectively mitigates this issue.The anti-windup species W 2 , shown in the inset figure, intervenes during the severe disturbance, preventing Z 1 from accumulating and, therefore, avoiding windup.Consequently, glucose levels never overshoot the healthy range.

DISCUSSION
The phenomenon of windup presents a substantial obstacle for any control system that aim to robustly regulate a target variable through integral feedback control.This challenge is especially evident in biomolecular and biomedical contexts where the limitations of actuators, sensors, and changes in processes may result in extended periods of poor performance or even loss of stability.With this in mind, our research focused on designing and evaluating biomolecular antiwindup strategies capable of mitigating windup effects and enhancing the robustness and reliability of biomolecular control systems.
In this study, we introduced three anti-windup topologies specifically adapted for biomolecular control systems that harness the RPAachieving properties of the AIF controller and its PID extensions.These topologies were carefully designed to account for biomolecular constraints such as positivity, promoter saturation, or resource limitations.We initially introduced the topologies using a high-level phenomenological representation involving switches.Subsequently, a sequestration reaction-based CRN was introduced to realize the behavior of the switches.Sequestration motifs are a common fixture in biomolecular designs to realize numerous functionalities including subtraction modules ( 62), biomolecular input-output linear systems (63), biomolecular integral control (6), biomolecular derivative operators (64,65,66), PID controllers as CRNs (22,67,68,69), band-pass filters (70), achieving ultrasensitive response (27), nonlinear activation functions for biomolecular neural networks (71), and stabilization near unstable equilibria (72), among others.This is predominantly attributed to the capability of sequestration reactions to perform comparison operations.In our work, we harnessed this characteristic to realize biomolecular antiwindup strategies, further underlining the essential role of sequestration motifs in biomolecular circuit design.
By leveraging Theorems 1 and 2, we established a connection between these biomolecular topologies and classical control-theoretic concepts.This association was made by calculating transfer functions and block diagrams that are exact in the limit of strong sequestration rates.This approach differs from previous approximation methods that relied on linearizations around the fixed point (22,23,66,69).Specifically, in this study, we exploit timescale separation which has previously been used to derive reduced order models describing the dynamics of biomolecular controllers (20,(73)(74)(75)(76).It is noteworthy to point out that in (73), a theorem similar to Theorem 1 is presented.Samaniego et al. (73) establish that if there exists some time T > 0 before which the concentrations of the sequestration pair do not cross, then simple sequestration networks reduce to subtraction operations for all t ≥ T. In contrast, Theorem 1 establishes a different model reduction result that is valid for more general networks involving a fast sequestration reaction and under less restrictive conditions.The model reduction results in Theorems 1 and 2 are valid for the dynamics over any compact time interval and for the fixed points, and they require no assumptions except Lipschitz continuity.In a sense, the two approaches complement each other since (73) addresses t ≥ T assuming that there exists a T that satisfies certain conditions; whereas Theorem 1 addresses t ∈ [0, T] for any T ≥ 0.
The promising potential of our methods to mitigate windup was illustrated both theoretically and through numerical simulations.Our (A) to (C) illustrate the genetic realizations for topologies i, ii, and iii, respectively.each realization is achieved using two genes, one driven by a constitutive promoter and another by an inducible promoter.the proteins expressed by these genes, V 1 and V 2 , engage in an intein-splicing reaction and are subsequently degraded.the function of V 2 varies across the topologies, acting as an activator in (A), a protease in (B), and a repressor in (c).the strength of the constitutive promoter or plasmid copy number can be manipulated to adjust the anti-windup threshold.Ad, activation domain; dBd, dnA binding domain; intc/intn, intein c/n; deg, degradation tag. the protein domains are color-coded to match the genes from which they are expressed.(47,48,61) as the controlled process.Windup effects are tested under a severe transient disturbance in the eGP rate.the gray response represents nondiabetic individuals not controlled by our circuits, while the blue and red responses represent patients with diabetes regulated by the AiF controller with and without anti-windup circuitry, respectively.With anti-windup circuitry, glucose levels return smoothly to the healthy range.Without it, Z 1 accumulation inhibits Z 2 and insulin, causing prolonged high glucose levels.Parameters are primarily from (10), except for the scaling of μ and θ by 100 to increase the controller's speed. in addition, η v = η w = 10 4 μmol −1 hour −1 , and v 0 = 100w 0 = 0.01 μmol/hour.the inset plots across all panels show the dynamics of the anti-windup species which become active when thresholds are exceeded.We have δ v1 = δ v2 = δ w1 = δ w2 = 1, and h 1 , h 2 , g 1 , g 2 from Fig. 8 are all identity functions.
simulations spanned three increasingly complex processes, thereby demonstrating the versatility and robustness of the proposed antiwindup topologies.From a basic gene expression process modeling transcription saturation to a highly intricate whole-body model of a patient with type I which is FDA approved, our topologies succeeded in preventing windup and maintaining effective control.We also put forth genetic designs capable of implementing the anti-windup strategies, based on inteins which were previously used to successfully build integral controllers.Thanks to their flexibility and small size, inteins hold considerable potential for use in various biomolecular circuits that involve sequestration reactions, such as those in our anti-windup designs.One primary limitation of the introduced biomolecular antiwindup strategies involves the potential saturation of the antiwindup circuitry itself.This issue can arise from saturation of relevant promoters or constraints due to limited shared resources.As demonstrated in the numerical simulations (particularly the inset figures in Fig. 10), the anti-windup species remain inactive as long as the controller operates within the preset thresholds of normal integral mode.Only when these thresholds are surpassed do the anti-windup species become active to mitigate windup.As such, the overall system benefits from the favourable property of low burden imposed by the added anti-windup circuitry, as long as no windup occurs.It is worth mentioning that, when properly tuned, the burden imposed by the anti-windup circuitry is minimal but not zero, as the molecules implementing the mechanism are still produced even in the absence of windup.However, if not properly tuned, then this intervention itself may saturate, compromising the effectiveness of the augmented anti-windup circuitry.Nevertheless, even in such a case, the dynamic range of the integrator is expanded compared to a system without anti-windup circuitry, thereby still reducing the impact of windup.In scenarios where the imposed burden is excessive, it may be necessary to distribute the AIF controller and the anti-windup circuitry across multiple cellular populations, following the approach suggested in (29,30).It is important to note that limited shared resources may lead to saturations affecting various network components such as reference, sensing, actuation, and anti-windup reactions and may also couple the dynamics among these components.A simulation study is presented in figs.S2 and S3 to demonstrate that, even with additional burden imposed by resource sharing, the anti-windup reactions ensure well-behaved dynamics.It is also analytically shown in the section S6 that the symmetry of the sequestration motif used to realize the switches in the anti-windup circuitry provides two crucial properties: (i) It enhances the robustness of the threshold parameter against the influence of shared resources, and (ii) it facilitates tuning the threshold which is managed through a production reaction analogous to the setpoint tuning in the AIF controller [see (27) and Box 1 for a similar discussion on threshold tuning].These properties are particularly useful as the threshold is an important parameter that separates the normal mode from the anti-windup protection mode.
For the effective implementation of the anti-windup switches depicted in Fig. 6A via sequestration reactions, strong sequestrations are paramount.While milder sequestrations can still curtail windup, they introduce steady-state errors even under normal operations (refer to fig.S6).To address this, designers should ensure that sequestration reactions occur more rapidly than other reactions.This is expected to be observed with inteins, where intein-splicing reactions are considerably faster than gene expression processes (77).
Moving forward, the proposed anti-windup topologies represent a promising advancement toward more robust biomolecular control systems.The genetic implementation of these topologies could offer powerful tools for a range of applications, from gene regulation to personalized therapeutics.Moreover, the broader approach of designing safe control systems with consideration for the unique features and constraints of biomolecular processes holds potential to forge previously unidentified paths in the design of robust and effective biomolecular control systems.A potential avenue for future exploration could involve examining the impacts of windup phenomena and the corresponding anti-windup strategies in a stochastic context.

MATERIALS AND METHODS
All simulations were performed using MATLAB.The MATLAB code for generating these simulations is available at the following Github repository https://github.com/Maurice-Filo/Biomolecular-Antiwindup-Circuits.Deterministic dynamics were simulated using the ODE23s solver, while the stochastic simulations in the Supplementary Materials were conducted using the Gillespie algorithm (78).The simulation of the glucose response in Fig. 10C was carried out using the SimBiology toolbox in MATLAB.Proofs of the theorems and lemmas are provided in the Supplementary Materials.

Fig. 1 .
Fig. 1.Integral windup.(A)A graphical portrayal of integral windup in a biomolecular system.the closed-loop network comprises an arbitrary process in a feedback configuration with a controller network involving integral feedback control.the controller's objective is to drive the regulated output Y to a prescribed setpoint r at steady state, despite constant disturbances and uncertainties-a property referred to as RPA. the controller operates by sensing the output Y, computing the integral of the error e(t) ≜ r − y(t), and feeding the integral control action back into the process.At steady state, ̇z = 0 and thus the error converges to zero.however, saturations in the actuator and/or sensor might obstruct this process, preventing the error from reducing to zero.this results in the integral of the error accumulating, leading to windup.this causes certain molecular species, denoted as Z, to increase uncontrollably, symbolized by the spiral "winding up" due to bouncing off of a saturation wall.therefore, equipping the controller with anti-windup circuitry becomes indispensable.(B) A typical behavior in the absence of anti-windup circuits.tolerable disturbances Δ can be promptly rejected by the controller, ensuring RPA.however, when the disturbance exceeds the manageable threshold, actuation u hits its saturation limit, and a control species Z starts to accumulate, triggering windup.even when the disturbance subsides, the error takes an extended period to return to zero because z(t) is substantially beyond the safe zone, resulting in a prolonged "unwinding" process.(C) Anti-windup circuits offer a solution by stepping in only when necessary, preventing the controller from reaching a dangerous state.this intervention introduces a smaller steady-state error, helps avoid the windup process, and ensures that, once the severe disturbance subsides, the error can quickly return to zero.

Fig. 2 .
Fig. 2. Classical integral feedback control.(A)Process or open-loop system.the input and output of the process are respectively denoted by u and y, while Δ denotes a disturbance or uncertainty.the dynamics of the process are described by a (nonlinear) dynamic operator Δ, and its input/output steady-state map for a constant disturbance or uncertainty Δ is denoted by  Δ .(B) Gene expression as an example of a simple process.the process is composed of two molecular species X 1 and X 2 representing the mRnA and proteins whose degradation rates are denoted by γ 1 and γ 2 , respectively.note that the translation rate is denoted by k 1 , and Δ represents a disturbance entering the dynamics as a basal transcription rate which could be due to a mutation in the promoter or constitutive expression from a duplicate gene.the input to this open-loop system is an induced transcription rate denoted by u, while the output is the protein concentration y = x 2 .the dynamics and input/output steady-state map are given here.(C) closing the loop with the classical iFc.ideally, the integral controller senses the output y and computes the integral of the error e(t) ≜ r − y, which signifies the deviation of the output from the setpoint r, multiplies this integral by a gain K I , and feeds the result back to the process through the actuator.in more practical scenarios, however, sensor and actuator saturation may occur, impeding the faithful transmission of signals between the controller and the process.the equations describing the dynamics of the closed-loop system are presented in eq. 9.

Fig. 4 .
Fig. 4. AIF control.(A)An arbitrary process regulated by the AiF controller.the process comprises L species X 1 through X L , with X L being the regulated output.the AiF controller comprises two species, Z 1 and Z 2 , which sequester each other at a rate η.Sensing is realized via the catalytic production of species Z 2 by the output X L at a rate θx L , whereas actuation is realized via the production of the input species X 1 at a rate u = kz 1 .the setpoint is encoded in the constitutive production of Z 1 at a rate μ.(B) Block diagram describing the exact closed-loop dynamics in the limit as η → ∞.By defining the setpoint r and integral gain K I , and an intermediate variable v, this diagram is generated.the diagram demonstrates that, in the strong sequestration limit, the AiF controller is exactly the classical iFc depicted in Fig.2c, excluding sensor saturation but including a one-sided actuator saturation due to the inherent positivity of the AiF motif.(C) integral windup induced by the positivity of the AiF controller. in this simulation example, the process is considered to be the gene expression model presented in Fig.2B. the green and gray shaded areas represent the set of admissible setpoints  and disturbances  r , respectively (see example 1 for more details).As the simulation demonstrates, the AiF controller maintains RPA as long as disturbances remain within admissible limits.however, if the disturbance shift to an inadmissible level, then the AiF controller falls short, exhibiting poor dynamic performance even after the disturbance returns to an admissible level-a clear demonstration of integral windup.numerical values: μ = 10, η = 100, k = θ = γ 1 = γ 2 = k 1 = 1 and Δ is varied.

Fig. 5 .
Fig. 5. AIF control with actuator and sensor saturation.(A) closed-loop network.the network structure remains the same as in Fig.4A, except the actuation and sensing rates are now represented as u = kh a (z 1 ) and θh s (y), respectively, to model nonlinear, monotonically increasing, and saturating functions.the input setpoint is designated as r in ≜ μ/θ, while the output setpoint, r out , is determined by the sensing function h s as per eq.19.(B) Block diagram describing the exact closed-loop dynamics in the limit as η → ∞. this diagram is similar to that in Fig.4B, with the actuator and sensor saturation modules replaced by the functions ψ a and ψ s , respectively.(C and D) integral windup induced by actuator/sensor saturation.the gene expression model from Fig.2Bis used here as the process.(c) tests setpoint tracking capabilities when the output setpoint, adjusted via μ, becomes transiently inadmissible.in contrast, (d) evaluates disturbance rejection properties by maintaining a fixed setpoint while varying the disturbance Δ until it becomes transiently inadmissible.the green shaded areas denote the set of admissible output setpoints  out . the gray shaded areas represent the set of admissible input setpoints (scaled by θ) and disturbances in (c) and (d), respectively (refer to example 3 for further details).the simulations reveal that the AiF controller maintains RPA effectively as long as the setpoint and disturbance stay within admissible limits.however, when either transitions to inadmissible levels, windup emerges manifesting as substantial accumulation of either z 1 or z 2 .numerical values: k 1 = γ 1 = γ 2 = 1, η = 100, θ = 15, κ a = κ s = 5, k = 8. in (c), Δ = 5 is fixed while μ is varied; whereas in (d), μ = 10 is fixed while Δ is varied.

Fig. 6 .
Fig.6.Anti-windup strategies.(A) three anti-windup topologies.the leftmost diagram depicts a closed-loop network regulated by one of three anti-windup topologies. in all topologies, the controller operates normally as an AiF controller ensuring RPA, as long as Z 1 and Z 2 levels stay below a threshold.Otherwise, the anti-windup mechanisms (red and blue switches) activate, shifting the controller to protection mode.topology i increases production of Z 1 in response to high Z 2 and vice versa.topology ii degrades Z 1 and Z 2 when their levels are high.topology iii reduces production rates of Z 1 and Z 2 when their concentrations exceed a threshold.the switches are represented by monotonic functions h 1 and h 2 .For topologies i and ii, these functions are zero until the threshold and then increase.For topology iii, they are one until the threshold and then decrease.Mathematical examples of these functions, realizable as cRns, are provided in the lower gray box (see Fig.7).introducing the intermediate variable v, input setpoint r in , integral gain K I , and functions ψ a and ψ s , as shown in the top gray box, generates the block diagram in the bottom left, depicting the reduced dynamics similar to Fig.5Bwith added anti-windup dynamics.(B) compact block diagram for topologies i and ii.these topologies use a conditional integration approach, where the controller integrates when v is within the safe range [−kβ i , kβ j ] and transitions to a leaky integrator outside this range to counter windup.topology i is associated with (i, j) = (1,2), while topology ii with (2,1).(C) Block diagram for topology iii, implementing reference/sensor conditioning.this reduces the input reference or sensor signal when v is outside the safe range [−kβ 2 , kβ 1 ].Specifically, high positive v scales down the input setpoint r in , while high negative v scales down the sensed signal w.

Fig. 8 .
Fig.8.CRN realizations of three anti-windup topologies.(A) topology i, (B) topology ii, and (C) topology iii. in each topology, sequestration networks, which consist of species (V 1 , V 2 ) and (W 1 , W 2 ) that sequester each other, are used to implement switch-like behaviors that intervene to prevent windup when necessary.topology i uses activator switches and so h 1 and h 2 are monotonically increasing beyond the threshold.topology ii uses a similar mechanism, but the switches catalyze the selfdegradations of Z 1 and Z 2 .here, h 1 (v 2 ) and h 2 (w 2 ) represent degradation rates that are monotonically increasing beyond the thresholds.the corresponding propensity functions, however, are given by h 1 (v 2 )z 1 and h 2 (w 2 )z 2 , respectively.As a result, if h 1 and h 2 are linear functions beyond the threshold, then the corresponding propensities are quadratic.last, topology iii introduces repressor switches, and so, h 1 (v 2 ) and h 2 (w 2 ) monotonically decrease beyond the thresholds.note that while the functions h 1 , h 2 , g 1 , and g 2 are kept intentionally general in this depiction, specific examples are given in (B) and (c). in (B), the functions are presented in their linear form, while in panel (c), h 1 takes on a hill form.

Fig. 9 .
Fig. 9. Intein-based genetic realizations of anti-windup topologies.(A) to (C) illustrate the genetic realizations for topologies i, ii, and iii, respectively.each realization is achieved using two genes, one driven by a constitutive promoter and another by an inducible promoter.the proteins expressed by these genes, V 1 and V 2 , engage in an intein-splicing reaction and are subsequently degraded.the function of V 2 varies across the topologies, acting as an activator in (A), a protease in (B), and a repressor in (c).the strength of the constitutive promoter or plasmid copy number can be manipulated to adjust the anti-windup threshold.Ad, activation domain; dBd, dnA binding domain; intc/intn, intein c/n; deg, degradation tag. the protein domains are color-coded to match the genes from which they are expressed.

Fig. 10 .
Fig.10.Simulations of anti-windup circuitry with three processes of increasing complexity.(A) the first simulation revisits the gene expression process, implementing the cRn realization of topology i to mitigate windup effects observed in Fig.5 (c and d).Anti-windup parameters are η v = η w = 100, v 0 = 10, w 0 = 20.(B) the second simulation illustrates compatibility with biomolecular Pid controllers [third-order Pid(22)] with saturations to regulate a complex process of six species.the plots show improved dynamic performance of the output X 6 .Parameter values from (22) are used, with additional parameters introducing saturations, basal expressions, and antiwindup circuitry: μ = 6, k = 15, θ = 300, κ a = κ s = κ p = 1000, b = 0.2, v 0 = w 0 = 3500, and η v = η w = η = 10 4 .(C) the final simulation uses a whole-body model of patients with type i diabetes(47,48,61) as the controlled process.Windup effects are tested under a severe transient disturbance in the eGP rate.the gray response represents nondiabetic individuals not controlled by our circuits, while the blue and red responses represent patients with diabetes regulated by the AiF controller with and without anti-windup circuitry, respectively.With anti-windup circuitry, glucose levels return smoothly to the healthy range.Without it, Z 1 accumulation inhibits Z 2 and insulin, causing prolonged high glucose levels.Parameters are primarily from(10), except for the scaling of μ and θ by 100 to increase the controller's speed. in addition, η v = η w = 10 4 μmol −1 hour −1 , and v 0 = 100w 0 = 0.01 μmol/hour.the inset plots across all panels show the dynamics of the anti-windup species which become active when thresholds are exceeded.We have δ v1 = δ v2 = δ w1 = δ w2 = 1, and h 1 , h 2 , g 1 , g 2 from Fig.8are all identity functions.